Submitted to ApJ 

Preprint typeset using I4TgX style emulateapj v. 10/09/06 



RAPIDLY ROTATING SUNS AND ACTIVE NESTS OF CONVECTION 

Benjamin P. Brown 

JILA and Dept. Astrophysical & Planetary Sciences, University of Colorado, Boulder, CO 80309-0440 

Matthew K. Browning 1 

Dept. of Astronomy, University of California, Berkeley, CA 94720-3411 

Allan Sacha Brun 

DSM/IRFU/SAp, CEA-Saclay & UMR AIM, CEA-CNRS-Universite Paris 7, 91191 Gif-sur-Yvette, France 

Mark S. Miesch 

High Altitude Observatory, NCAR, Boulder, CO 80307-3000 
AND 
JURI TOOMRE 

JILA and Dept. Astrophysical & Planetary Sciences, University of Colorado, Boulder, CO 80309-0440 

(Received April 30, 2008; Revised July 1, 2008; Accepted August 12, 2008) 

Submitted to ApJ 

ABSTRACT 

In the solar convection zone, rotation couples with intensely turbulent convection to drive a strong differential 
rotation and achieve complex magnetic dynamo action. Our sun must have rotated more rapidly in its past, as is 
suggested by observations of many rapidly rotating young solar-type stars. Here we explore the effects of more 
rapid rotation on the global-scale patterns of convection in such stars and the flows of differential rotation and 
meridional circulation which are self-consistently established. The convection in these systems is richly time 
dependent and in our most rapidly rotating suns a striking pattern of localized convection emerges. Convection 
near the equator in these systems is dominated by one or two nests in longitude of locally enhanced convection, 
with quiescent streaming flow in between at the highest rotation rates. These active nests of convection maintain 
a strong differential rotation despite their small size. The structure of differential rotation is similar in all of 
our more rapidly rotating suns, with fast equators and slower poles. We find that the total shear in differential 
rotation A57 grows with more rapid rotation while the relative shear Af2/f2o decreases. In contrast, at more 
rapid rotation the meridional circulations decrease in energy and peak velocities and break into multiple cells 
of circulation in both radius and latitude. 

Subject headings: convection - Sun:interior - Sun:rotation - stars:interiors - stars:rotation 



1 . CONVECTION, ROTATION AND MAGNETISM 

Our sun is a magnetic star whose cycles of magnetic activ- 
ity must arise from organized dynamo action in its interior. 
This dynamo action is achieved by turbulent plasma motions 
in the solar convection zone, which spans the outer 29% of the 
sun in radius. Here vigorous convective motions and rotation 
couple to drive the differential rotation and meridional circu- 
lation. These global-scale flows are important ingredients in 
stellar dynamo theory, providing shear which may build and 
organize fields on global scales. When our sun was younger, 
it must have rotated much more rapidly, as is suggested both 
by the solar wind which continually removes angular momen- 
tum from the sun and by many observations of rapidly rotat- 
ing solar-like stars. In more rapidly rotating suns the coupling 
between rotation and convection is strong and must continue 
to drive global scales of flow. Understanding the nature of 
convection, differential rotation and meridional circulation in 
more rapidly rotating stars is a crucial step towards under- 
standing stellar dynamos. 

The manner in which the sun achieves global-scale dynamo 
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action is gradually being sorted out. Helioseismology, which 
uses acoustic oscillations to probe the radial structure of the 
star as well as the convective flows beneath the surface, has 
revealed that the solar differential rotation profile observed at 
the surface prints throughout the convection zone with two 
prominent regions of radial shear. The near-surface shear 
layer occupies the outer 5% of the sun, whereas a tachocline 
of shear at the base of the convection zone separates the strong 
differential rotation of that zone fro m the uniform rotatio n 
of the deeper radiative interior (e.g., [Thompson et alj|2003l) . 
The solar global magnetic dynamo, responsible for the 22- 
year cycles of activity, is now believed to be seated in that 
tachocline. In such interface dynamo models, magnetic fields 
generated in the bulk of the convection zone are pumped into 
the tachocline where the radial shear builds strong toroidal 
magnetic fields, with magnetic buoyancy leading to loops 
that rise upward and erupt throu gh the solar surface (e.g., 
ICharbonneaull2005l : lMieschll2005l) . The differential rotation 
plays an important role in building and organizing the global- 
scale fields while the meridional circulations may be impor- 
tant for returning flux to the base of the convection zone and 
advecting it equatorward, enabling cycles of magnetic activ- 
ity. By studying the coupling of rotation and convection over a 
range of conditions, we are likely to learn about the operation 
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of the current solar dynamo and about the nature of dynamos 
operating in our sun's past and in other solar-like stars. 

Observations of young solar-like stars indicate that they 
rotate as much as 50 times faster than the current solar 
rate. Many of these more rapidly rotating suns possess 
strong magnetic fields. A correlation between rotation and 
magnetic activity is observed over a range of stellar types 
and populations, indicating that more ra pidly rotating stars 
may have stronger st ellar dynamos (e.g., iNoves et al.lll984t 
Patten & Simon jl996l) . Probing the nature of these dynamos 
and the impact of faster rotation on the internal stellar dy- 
namics requires both accurate observations and detailed dy- 
namical models of the stellar interiors. The faster flows of 
differential rotation are much easier to detect than the rela- 
tively slow motions associated with meridional circulations; 
observations across the HR diagram indicate that differential 
rotation is a common feature in many stars. Asteroseismic ob- 
servations with the Kepler and Corot missions may soon be- 
gin to constrain the internal rotation structure. At present only 
measurements of surface differential rotation are available, as 
assessed w ith a variety of techniques including p hotometric 
variabilit y (|Donahue et1d1ll996l:IWalker et alj2007l) . Doppler 
imagin g (Donat i et alJl2003l) and Fourier transform methods 
dReiners & Schmittll2003l ). 

Advances in supercomputing have enabled 3-D simulations 
that are beginning to capture many of the dynamical ele- 
ments of the solar convection z one. Ear l y global-sca l e sim- 
ulations of solar convection by iGilmanl (1 19751 1 19771 1 1979b 
under the Boussi nesq approximation were exte nded by the pi- 
oneering work of Gilman & Glatz maierl (1 1 98 lh . Such global- 
scale simulations of solar convection conducted in full spher- 
ical shells sought to capture the largest scales of convective 
flows and began to study how they can establish differential 
rotation and meridional circulations. However, the range of 
spatial and temporal scales present in solar convection are 
vast and thus the computational resources required by the 
modeling are daunting. Through recent advances in mas- 
sively parallel computer architectures, solar convection sim- 
ulations are now beginning to make detailed contact with the 
observational constraints provided by helioseism ology (e.g., 
iBrun & Toomrel 120021: iMiesch et al.ll2006l 120081) . Other ef- 
forts have focused on the vigorous turbulence and the dy- 
namo action achie ved in the bulk of the solar convection zone 
(IBrun et al.ll2004l) . with recent studies beginning to include 
the tachocline as a region of p enetrative overshoot, s hear, and 
magnetic field amplification (Bro wning et al.ll2006l) . Facili- 
tated by these computational advances, models of convection 
and dynamo action within the cores of A-type sta r s have also 
begun to be investigated dBrowning et alj l2004t IBrun et all 
120051: iFeafherstone et alJl2007h . 

To date, most models of stellar differential rotation in 
stars like our sun that rotate more rapidly have been car- 
ried out in 2-D un der the simplifying assumptions of mean 
field theory (e.g iRudjre reLaLl 119981: iKuker & Stixl 1200 It 
iKiiker & Riidi gerl 12005 alibi) . The time is ripe to pursue the 
question with fully 3-D simulations of global-scale stellar 
convection. 

In order to study solar-like stars that rotate faster, our pre- 
vious work focused on a series of 3-D compressible simu- 
lations within a spherical shell using the anelastic spherical 
harmonic (ASH) code for stars rotating f rom one to five times 
the current solar rate (Brown et al. 2004). These preliminary 
hydrodynamic simulations explored how stellar convection 
changes with more rapid rotation, including the differential 



rotation and meridional circulation that is achieved. Com- 
parable studies have been carried out by iBallot etafl (120071) 
studying younger stars with deeper convection zones. In our 
simulations we found that remarkable nests of vigorous con- 
vection emerge in the equatorial regions. Namely, convective 
structures at low latitudes about the equator can exhibit strong 
spatial modulation with longitude, and at high rotation rates 
the convection is confined to narrow intervals (or nests) in 
longitude. In the more turbulent simulations presented in this 
paper, which span a larger range of rotation rates from one to 
ten times the current solar rate, the phenomena of modulated 
convection has persisted and the active nests of convection are 
prominent features at the higher rotation rates. These nests of 
localized convection persist for long intervals of time and de- 
spite their small filling factor maintain a strong differential 
rotation. 

The emergence of spatially localized convective states has 
been observed in other systems, particularly in theoretical 
studies of doub ly-diffusive systems such as therm osolutal 
convection (e.g. JSpina et al.lfl9 98: Batiste et al. 2006), in lab- 
orator y studies of convection in binary fluids (e.g. JSurko et alj 
1199 lh . and in simulations of magnetoc onvection where iso - 
lated "convectons" have been observed ( Blanchfl oweril 1 999h . 
In shells of rapidly rotating fluid, temporally intermittent 
patches of localized conve ction emerged in Bou ssinesq simu- 
lations of the geodynamo dGrote & Bussdl2000l) . In many of 
these systems, spatial modulation occurs in the weakly non- 
linear regime close to the onset of convection. In contrast, 
our simulations of solar convection are in a regime of fully 
developed turbulent convection. 

We describe briefly in $2] the 3-D anelastic spherical shell 
model and the parameter space explored by our simulations. 
In $3] we discuss the nature of convection realized in more 
rapidly rotating stars and the emergence of spatially-localized 
patterns of convection. In §^4]|6]we examine the global-scale 
flows realized in our simulations, including differential rota- 
tion and meridional circulation, and their scaling with more 
rapid rotation. A detailed exploration of the active nests of 
convection is presented in $71 We reflect in §|8]on the signifi- 
cance of our findings. 

2. TOOLS TO MODEL 3-D STELLAR CONVECTION 

To study through simulations the global scales of motion 
likely achieved in stellar convection zones, we must employ 
a global model which simultaneously captures the spherical 
shell geometry and admits the possibility of zonal jets and 
large eddy vortices, and of convective plumes that may span 
the depth of the convection zone. The solar convection zone 
is intensely turbulent and molecular values of viscosity and 
thermal diffusivity in the sun are estimated to be very small. 
As a consequence, numerical simulations cannot hope to re- 
solve all scales of motion present in real solar convection and 
a compromise must be struck between resolving the important 
dynamics within small regions and capturing the connectivity 
and geometry of the global scales. We opt here for the latter 
by studying a full shell of convection. 

2.1. Anelastic Simulation Approach 

Our tool for exploring stellar convection is the anelastic 
spheri cal harmonic (ASH ) code, which is describe d in de- 
tail in IClune et alJ (1 19991) and in IBrun et all (120041) . ASH 
is designed to run efficiently on massively parallel architec- 
tures. ASH solves the three-dimensional (3-D) anelastic equa- 
tions of motion in a rotating spherical shell using the pseudo- 
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spectral method. The thermodynamic variations and the three 
components of mass flux are expanded in spherical harmonics 
to resolve their horizontal structure and in Chebyshev polyno- 
mials to resolve their radial structure. The anelastic approx- 
imation is used to capture the effects of density stratification 
without having to resolve sound waves which have short pe- 
riods (about 5 minutes) relative to the dynamical time scales 
of most interest (months to years). The mass flux remains 
divergence-free by using a poloidal-toroidal representation. 
Temporal discretization is accomplished using a semi-implicit 
Crank-Nicholson time-stepping scheme for linear terms and 
an explicit Adams-Bashforth scheme for nonlinear terms. Un- 
der the anelastic approximation the thermodynamic variables 
are linearized about their spherically symmetric and evolving 
mean state with density p, pressure P, temperature f and spe- 
cific entropy S. Fluctuations about this state are denoted as 
p, P, T and S respectively. In the uniformly rotating refer- 
ence frame of the star, the resulting hydrodynamic equations 
expressing the conservation of mass, momentum and internal 
energy are in turn: 

V-(pv) = 0, (1) 



dv -> 
— + (v- V)v+2S1 () x v 
at 

-V(P+P) + (p+p)g-W-V, 



-dS 



pT—= V- \K r pc p V(T+T) + K pTVS+KpTVS 
at 1 J 



-pTv- V{S+S) + 2pv 



(2) 



(3) 



where v = (v r , ve , V^) is the local velocity in the rotating frame 
of constant angular velocity Oq> S is tne gravitational accel- 
eration, c p is the specific heat at constant pressure, n r is the 
radiative diffusivity and D is the viscous stress tensor, given 
by 

V ij = -2pv e y -i(V-v)*y » (4) 

where e,j is the strain rate tensor. Here v and k are the dif- 
fusivities for vorticity and entropy. We assume an ideal gas 
law 

P=npT, (5) 

where 1Z is the gas constant, and close this set of equations 
using the linearized relations for the thermodynamic fluctua- 
tions 

p _P T _ P S 
~p Z ~P~T~rP~V p 

The mean state variables are evolved with the fluctuations, 
thus allowing the convection to modify the entropy gradients 
which drive it. 

ASH is a large-eddy simulation (LES) code, with subgrid- 
scale (SGS) treatments for scales of motion which fall be- 
low the spatial resolution in our simulations. We treat these 
scales with effective eddy diffusivities, v and n, which repre- 
sent momentum and heat transport by unresolved motions in 
the simulations. Here for simplicity v and k are taken as func- 
tions of radius alone, and proportional to p" 1 / 2 . These simula- 
tions a re similar to case AB as reported in Bru n & Toomrd 
(2002) though with a different SGS functional form (there 
v,k oc p~ l ), and we shal l consider faster ffy. Our adopted 
SGS v ariat ion here, as i nlBrun et al.l (120041) . [Browning et al.l 
(2006) and Ball ot et al.l (120071) . yields lower diffusivities near 
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the top of the layer and thus higher Reynolds numbers. The 
thermal diffusion acting on the mean entropy gradient kq is 
treated separately and occupies a narrow region in the upper 
convection zone. Its purpose is to transport heat through the 
outer surface where radial convective motions vanish. 

Our simulations are still separated by many orders of 
magnitude from the intensely turbulent conditions present 
within the solar convection zone. They are likely to cap- 
ture many aspects of the dynamics of solar convection, and 
we a re encouraged by the success that similar simulations 
(e.g. IMiesch et al.ll2000l iBrun & Toomrell2002t iMiesch et ail 
2006, 2008) have had in beginning to match the detailed 
observational constraints for differential rotation within the 
solar convection zon e provided by helioseismology (c.f. 
iThompson et~ai1l2003l) . 

We shall impose velocity boundary conditions which im- 
pose no net torques on the shell. These are in turn: 

1 . Impenetrable top and bottom: v r = 0, 

2. Stress-free top and bottom: 

(d/dr)(v 6 /r) = (d/dr)(v^/r) = 0, (7) 

3. Constant entropy gradient at top and bottom: dS/dr = 
const. 

Some of our simulations are initialized by perturbing a qui- 
escent state in solid body rotation. The growth of convection 
leads to velocity correlations that serve to redistribute angu- 
lar momentum within the shell, building a differential rotation 
and meridional circulation. We evolve the simulation for long 
periods compared variously to convective overturning times, 
rotation periods or typical diffusive times. Other simulations 
were started from these evolved states and then run for long 
epochs after all adjustments have been made to the frame ro- 
tation rate Qq and to the viscous and thermal diffusivities. 

2.2. Studies of Global-Scale Convection in G-type Stars 

Our numerical model is a relatively simple description of 
the solar convection zone that captures the essential spherical 
geometry and global connectivity of that domain. Solar val- 
ues are taken for heat flux, mass and radius and a perfect gas 
is assumed. Near the solar surface the H and He ionization 
zones, coupled with radiative losses, drive intense convec- 
tive motions on very small scales which appear at the surface 
as granulation. Capturing granulation in a global simulation 
would require spherical harmonic degrees of order 4000 and 
this is currently too demanding. We therefore position the top 
of our domain slightly below these ionization layers. In these 
simulations we also omit the stably stratified radiative interior 
and the shear layer at the base of the convection zone known 
as the tachocline. We focus here on the bulk of the convection 
zone, with our computational domain extending from 0.72/? Q 
to O.96/?0, thus spanning 172 Mm in radius. The reference or 
mean state of our thermodynamic varia bles is derived f rom a 
one-dimensional solar structure model (Brun et al. 2002]) and 
is continuously updated with the spherically symmetric com- 
ponents of the thermodynamic fluctuations as the simulations 
proceed. These values are illustrated in Figure Q] after convec- 
tion has readjusted the stratification. 

Our studies here explore a variety of solar-like stars rotat- 
ing from 1 to 10 £Iq (cases G1-G10). All cases use the same 
initial stellar structure. We seek here to explore the general 
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TABLE 1 

Parameters for Primary Simulations 



Case 




Ra 


Ta 


Re 


Re' 


Ro 


Roc 


;/ 


K 




Gl 


96 x 256 x 512 


3.22 X 10 4 


3.14 x 10 5 


84 


63 


0.92 


0.61 


2.75 


11.0 


1 


G2 


96 x 256 x 512 


1.75 X 10 5 


3.21 x 10 6 


205 


85 


0.55 


0.45 


1.72 


6.87 


2 


G3 


96x 256 x512 


4.22 X 10 5 


1.22 x 10 7 


326 


103 


0.41 


0.36 


1.32 


5.28 


3 


G4 


96x256x512 


7.89 X 10 s 


3.18 x 10 7 


433 


119 


0.33 


0.31 


1.09 


4.36 


4 


G5 


96x256x512 


1.29 X 10 6 


6.70 x 10 7 


543 


133 


0.28 


0.27 


0.94 


3.76 


5 


G7 


192 x 512 x 1024 


2.63 x 10 6 


2.06 x 10 8 


763 


154 


0.22 


0.22 


0.75 


3.01 


7 


G10 


192 x 512 x 1024 


5.58 x 10 6 


6.74 x 10 s 


1051 


188 


0.17 


0.18 


0.59 


2.37 


10 


G3a 


96x256x512 


7.83 x 10 5 


2.41 X 10 7 


528 


158 


0.50 


0.34 


0.94 


3.76 


3 


G3b 


192 x256x512 


2.26 x 10 6 


8.02 X 10 7 


1121 


324 


0.70 


0.32 


0.52 


2.06 


3 


G5b 


192 x 512 x 1024 


4.03 x 10 s 


2.23 X 10 8 


1347 


274 


0.41 


0.26 


0.52 


2.06 
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NOTE. — All simulations have inner radius rt, ot = 5.0 X 10 cm and outer radius of r top = 6.72 X 10 cm, with 
L = (/top - Ijoi) = 1-72 X 10 10 cm the thickness of the spherical shell. Evaluated at mid-depth are the Rayleigh 
number Ra = (-dp/dS)(dS/dr)gL* / pun, the Taylor number Ta = 4£2 2 ,L 4 /i/ 2 , the rms Reynolds number Re = 
VrmsL/u and fluctuating Reynolds number Re' = v' Ims L/u, the Rossby number Ro = u/2Qq , and the convective 
Rossby number Roc = (Ra/TaPr)' /2 . Here the fluctuating velocity v' has the differential rotation removed: v' = v- 
(v), with angle brackets denoting an average in longitude. The Prandtl number Pr = u/k is 0.25 for all simulations. 
The viscous and thermal diffusivity, v and k, are quoted at mid-depth (in units of 10 12 cm 2 s~'). The rotation rate 
of each reference frame Qq is in multiples of Qq = 2.6 X 10~ 6 rad s or 414 nHz. The viscous time scale at 
mid-depth t„ = L 2 /u is 1250 days for case Gl and is 3640 days for case G5. Additional cases considered other 
rotation rates at 1.25, 1.5, 1.75 and 6 Qq. 

convection scales with rotation as Ra c o c f^ 3 for Boussines q 
convection (e.g., Chandrasek fiarl fl96lt iDormv et alj [2004). 
Lower diffusivities lead to both longer viscous and thermal 
diffusion time scales and to flows possessing finer spatial 
scales. Achieving equilibrated states in these systems re- 
quires high resolution simulations carried out over extended 
periods. We have taken a middle ground between attempt- 
ing to maintain constant supercriticality (which may require 
scaling j/,k oc f^ 2 ) and keeping resolution requirements rea- 

-2/3 

sonable by scaling our diffusivities as v, k oc fl . However, 
all of our cases here are still highly supercritical, noting that 
the critical Rayl eigh number for these simulations at 1 is 
Ra f - 1000 dGilman & Glatzmaierlll981l:lMieschlll998l) . We 
have maintained a constant Prandtl number Pr = vj k = 0.25 in 
all of our simulations. The parameters of our models are de- 
tailed in TableQ] Our choice of scalings for the eddy diffusiv- 
ities with rotation rate may have some influence on the nature 
of the convective patterns and mean flows we achieve. To as- 
sess some of the sensitivity the choice of our path through pa- 
rameter space, we have also sampled a limited range of more 
turbulent simulations at a few rotation rates (cases G3a, G3b 
and G5b). 

All of the simulations discussed here are at approximately 
the same level of maturity in their evolution. Case Gl was the 
progenitor case at 1 f2 Q and was evolved for some 3000 days 
after b ranching away from case AB from Bru n & Toomrd 
(2002), which itself has seen about 10000 days of total simu- 
lated life. Starting with this case, each subsequent simulation 
was spun up from the next fastest case (i.e., G3 was spun up 
from G2) and evolved for over 4000 days, or many hundreds 
of rotation periods. At this point, all cases appear to be sta- 
tistically stationary in terms of the angular momentum fluxes 
and the kinetic energies. We believe that the differential rota- 
tion profiles presented are effectively stationary, though there 
are small fluctuations as determined from short averages over 
a few rotation periods. Certain cases (including G5) were 
evolved for much longer intervals (over 10000 days and more 
than 2000 rotation periods) to explore the long-term behavior 
of convective patterns in these rapidly rotating systems. To 




0.75 0.80 0.85 0.90 0.95 

FIG. 1 . — Radial variation of mean stellar structure in the ASH models, (a) 
Entropy gradient (dS/dr) for cases Gl, G5 and G10 (as labeled). At higher 
rotation rates the entropy gradient becomes steeper throughout the convection 
zone, even for our most turbulent cases (case G5b, long dashes), (b) Temper- 
ature and density (latter ranging from 0.203 to 0.008 gcnT 3 in the region 
simulated) for case Gl. (c) Pressure scale height Hp (in Mm and fractional 
solar radii), for case Gl, with cases G2-G10 similar in their mean stratifica- 
tion. 



effects of rotation on stellar convection rather than the evo- 
lution of a particular star, which would require modifications 
to the stellar structure as the star aged. In surveying the ef- 
fects of more rapid rotation on global-scale convection, we 
sought to achieve reasonably high levels of turbulence in the 
resulting flows. Thus our trajectory through the parameter 
space of £lo,^, and k attempts to maintain strong nonlinear- 
ity without having increasing the Hq laminarize the convec- 
tion. As we increased the rotation rate, we simultaneously 
decreased the effective eddy diffusivities v and k to main- 
tain the supercriticality of the simulated turbulent convection. 
We note that the critical Rayleigh number for the onset of 
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FIG. 2. — Changes in convective patterns with increasing rotation rate for 
mildly turbulent cases at (a) one, (b) three and (c) rive times the solar rate. 
Shown as snapshots are radial velocities near top of layer in global Mollweide 
projection, with upflows light and downflows dark. Poles are at top and bot- 
tom, and the equatorial region appears at middle, with equator indicated by 
bold dashed line. Thin dashed lines denote circles of constant latitude or 
longitude, and the thin surrounding line indicates the location of the stellar 
surface at Rq . A striking pattern of convection localized into nests near the 
equator emerges as the rotation rate increases. 

test that our results are not unduly subject to hysteresis in the 
system, we explored a branch of cases which were succes- 
sively spun down from 5 £Iq to 1 f2 Q . No significant hystere- 
sis was found. 

3. CONVECTION IN RAPIDLY ROTATING SUNS 
3.1. Early Results of Modulated Convection 

In our previous simulations of rapidly rotating suns we 
found that strongly local ized states of conv ection emerged 
with more rapid rotation (Bro wn et al.l l2004). A selection of 
these simulations in Figure [2] present snapshots of radial ve- 
locity near the top of the domain in global Mollweide projec- 
tion, showing the entire spherical surface with minimal dis- 
tortion. With more rapid rotation, a prominent longitudinal 
modulation appears in the patterns of equatorial convection. 
At the higher rotation rates the equatorial convection is con- 
fined to one or two nests, with streaming zonal flow filling the 
regions outside these nests of convection. These nests can per- 
sist for intervals spanning many hundreds of rotation periods, 
often with little change. Two nest states sometimes evolve 
into single nest states as one nest overtakes another. 

The simulations shown in Figure |2] are less turbulent than 
the cases presented in the rest of this paper, each possessing 
Reynolds numbers that are about three-fold smaller near the 
surface than in our new simulations. In the former models, 
a large portion of the energy transport in the upper convec- 
tion zone was carried by unresolved scales of motion. This 



parametrized SGS flux, represented by kq, dominated trans- 
port in the upper 30% of the convection zone, leading to 
weaker enthalpy transport and weaker resolved convection 
there. This parametrized flux acts as volume cooling term that 
removes flux from the regions where it is dominant; the dy- 
namics were influenced by the presence of this cooling layer. 
The cases presented in this paper have a narrower unresolved 
flux layer, confined now to the upper 10% of the convection 
zone, and consequently much more vigorous convection is re- 
alized throughout the domain. In these more turbulent cases, 
the phenomena of localized nests of convection is realized at 
somewhat higher rotation rates. 

3.2. Convective Patterns and Evolution with Rotation 

The variation of convective patterns with increasing rota- 
tion rate VLq in our more rapidly rotating suns is illustrated 
in Figure [3] Snapshots of the radial velocity near the top of 
the domain (0.95/?q) are shown in Mollweide projection for 
four cases: Gl, G3, G5 and G10. The convection patterns 
are complex and time dependent, with asymmetries between 
the upflows and downflows owing to the density stratification: 
narrow, fast downflow lanes are surrounded by broad, rela- 
tively weak upflows. 

There is a clear difference in both the scale and structure 
of convection at high and low latitudes. In the equatorial re- 
gions (roughly ±30° in latitude), the downflows organize into 
large structures (loosely called banana cells) aligned with the 
rotation axis, thus extending in the north-south direction. At 
high rotation rates this tendency for alignment becomes pro- 
nounced, largely in the spirit of the Taylor-Proudman theo- 
rem, and the downflow network exhibits little of the east-west 
branching visible in case Gl. These downflow lanes propa- 
gate in a prograde sense relative to the bulk rotation rate and 
do so more rapidly than the differential rotation which they 
themselves establish. The nests of convection, when they ap- 
pear at the higher rotation rates, propagate at an intermediate 
rate as denoted by the heavy dashed contours in Figure ^f-h. 
We defer discussion of the nature of these nests of convection 
until later. Individual convective cells persist for about 10 to 
30 days. 

In the higher latitude regions, the convection cells are more 
isotropic and the downflow network organizes on smaller 
scales. Convection in these regions is vigorous and com- 
plex, with upflows and their downflow networks in a constant 
dance. The convective cells have a cusped appearance, with 
downflows leading upflows as both propagate in a retrograde 
fashion (most apparent in Figs [3j>, c). Strong vortical plumes 
form in the interstices of the downflow network at both high 
and low latitudes. In the polar regions above the middle of 
the convection zone, the sense of vorticity in these downflow 
plumes is generally cyclonic: counterclockwise in the north- 
ern hemisphere and clockwise in the southern. As they de- 
scend through the mid-layer their vorticity changes and they 
become largely anti-cyclonic. In contrast, the polar upflows 
are anti-cyclonic at all depths. 

The latitudinal variation of convection patterns can be in 
part understood by considering a cylinder tangent to the 
base of the convection zone and aligned with the rotation 
axis. Within our geometry, this tangent cylinder intersects 
the outer boundary at about ±42° of latitude. It is well 
known that in a rotating convective shell the flow dynam- 
ics are different inside and outside of the tangent cylinder, 
owing to differences in the connectivity of the flows, the in- 
fluence of the Coriolis forces and distance from the rotation 



6 



Brown, Browning, Brun, Miesch & Toomre 




FIG. 3. — Radial velocity patterns in Mollweide projection at 0.95Rq (left) and differential rotation profiles (middle, right) with increasing rotation rate in (a,e) 
for case Gl, (b,f) for G3, (c,g) for G5, and (d,h) for G10. At higher rotation rates the horizontal scale of convective cells shrinks at all latitudes and cells are more 
strongly aligned with the rotation axis. A striking pattern of modulated convection emerges at low latitudes with faster rotation, consisting of spatially modulated 
or patchy convection. These active nests of convection are propagating structures which persist for long periods of time. At middle are profiles of mean angular 
velocity Q with radius and latitude. These differential rotation profiles all involve fast equators (prograde relative to the frame rate Qq, indicated by tickmark 
on scale) and a monotonic decrease of Q as the poles are approached. At right are radial cuts of the angular velocity at selected latitudes, as labeled. The dark 
dashed contour denotes the constant propagation rate of the nests where discernible. 
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axis (e.g., Buss^ ll970l) . These differences become more evi- 
dent as the rotation rate, and hence the rotational constraints 
on the convection, increases. With more rapid rotation the 
longitudinal extent of the convective cells becomes progres- 
sively smaller. Linear theory, in the Boussinesq approxi- 
mation, predicts that the wavenumber of the most unstable 

mode scales with rotation as m oc Ta 1 / 6 oc (e.g., 
lChandrasekhaJll96Ulbormv et alj|2004l) for polar and equa- 
torial convection. This effect is fo und in anelastic systems as 
well (Glatzmaier & Gilman|[l981h and in our simulations is 
evident at both high and low latitudes. 

Shown at right in Figures [3^-/1 are the profiles of differ- 
ential rotation (as angular velocity Jl) realized in these sim- 
ulations. These il(r, 8) profiles are azimuthally and tempo- 
rally averaged over a period of roughly 200 days. All of our 
more rapidly rotating solar-like stars exhibit solar-like dif- 
ferential rotation profiles, with prograde (fast) equators and 
retrograde (slow) poles. Contours of constant angular veloc- 
ity are aligned nearly on cylinders, influenced by the Taylor- 
Proudman theorem, though recent simulations of solar con- 
vection suggest that this is sensitive to the treatment o f the 
bottom thermal boundary condition (Mi esch et alJl2006l) . As 
first shown in mean-field models by Rempel (2005) an d the n 
in models of global-scale convection by Mie sch et al.l (2006), 
introducing a weak latitudinal gradient of entropy at the base 
of the convection zone, consistent with a thermal wind bal- 
ance in a tachocline of shear, can serve to rotate the Jl con- 
tours toward the more radial alignment deduced from helio- 
seismology without significantly changing either the overall 
£1 contrast with latitude or the convective patter ns. We ex- 
pect si milar behavior here, as briefly explored by Ball ot et all 
(120071) . but have not explored this issue in detail at the higher 
rotation rates. More rapidly rotating suns may very well also 
possess tachoclines, but at this stage there is no observational 
evidence of this. Thus we have simplified these simulations 
by imposing a constant entropy at the bottom boundary. Our 
contours of £1 in Figure[3]show some differences between the 
northern and southern hemispheres, particularly at higher lati- 
tudes, and these differences decrease with more rapid rotation. 
The patterns of convection are not simply symmetric about 
the equator, and thus the accompanying mean zonal flows can 
be expected to show some variations between the two hemi- 
spheres. Also shown are radial cuts of £1 at six fixed latitudes 
that make evident the angular velocity contrasts with radius 
and latitude achieved in these simulations. The absolute con- 
trast in latitude and radius clearly grows with rotation rate, 
and will be discussed in ^4] 

A most striking result of our simulations is the emergence 
of persistent, spatially modulated convection in the equatorial 
regions at high rotation rates. At these low latitudes, convec- 
tion becomes modulated in longitude and forms distinct active 
nests where the convective vigor is enhanced compared to the 
regions outside. The amplitude of the convective motions and 
enthalpy transport is larger within these nests, and indeed at 
the highest rotation rates, convection in the equatorial band is 
confined entirely to the nests. These nests of convection prop- 
agate at a velocity distinct from either the zonal flow of differ- 
ential rotation or that of the individual cells of convection and 
persist for very long periods of time (more than 5000 days in 
case G5). The nature of these active nests of spatially local- 
ized convection will be explored in detail in ^7] 

Weak modulation in longitude is already evident at low ro- 
tation rates. When long time series are considered we have 




FIG. 4. — Variation of mean zonal velocity (v,*) with latitude for case G5 
sampled at six radial cuts as labeled and shown here relative to the uniform 
propagation rate of the nests of convection. The nests experience a strong 
prograde zonal flow (positive) near the top of layer and a prominent retro- 
grade flow within the lower half. 

positively identified nests of convection in all simulations ro- 
tating at Oo > 3 0,0, As the rotation rate increases, the mod- 
ulation level gradually increases; at the highest rotation rates 
(> 7 fi Q ) the equatorial convection is almost solely confined 
to the nests. The convection realized in case G10 (Fig. [3]f) 
is marked by this extreme modulation, with strong upflows 
and downflows inside the nest and very little convection in 
the surrounding regions. These most rapidly rotating cases 
maintain a strong differential rotation profile, even though the 
equatorial convection occupies only a narrow interval in lon- 
gitude. The regions outside the nest are filled with fast stream- 
ing zonal flows consistent with the differential rotation. 

3.3. Radial Connectivity of Convection 

The nests of enhanced convection span the convection zone 
and propagate everywhere at a constant prograde angular ve- 
locity relative to the bulk rotation rate of the star. A contour 
corresponding to this characteristic propagation rate is over- 
plotted on the differential rotation profiles in Figures for 
cases G3, G5 and G10. As is evident from these profiles, 
the angular velocity associated with the differential rotation 
exceeds that of the propagation rate of the nests near the sur- 
face and is slower than that near the base of the convection 
zone. The nests of convection therefore live within an envi- 
ronment of substantial zonal shear with radius, as is quantified 
for case G5 in Figure |4] Here the shearing zonal velocity of 
differential rotation is plotted in latitude at six radial depths. 
At all depths there is substantial zonal flow through the nests 
of convection. 

In other studies of solar convection, we have found that 
strong downflo w lanes extend throughout the entire depth 
of the domain dMiesch et all 120001: ferun & Toomrel 120021: 
MieschetaL|E008). In our more rapidly rotating stars, this 
connectivity with depth changes markedly, as is illustrated for 
case G5 in Figure [5] showing radial velocities throughout the 
convection zone. In these rapidly rotating suns, the strong 
variation of mean zonal flow with radius in the equatorial re- 
gions prevents all but the strongest downflows from spanning 
the convection zone. Within the nest of enhanced convection 
the plumes are able to traverse the convection zone. Yet in the 
quieter regions outside, the weaker downflow plumes are trun- 
cated by shear before reaching the middle of the convection 
zone. It is evident (Fig |5j;) that the amplitude of convective 
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LongHude 

FIG. 5. — Connectivity of radial velocity with depth in case G5, shown at 
same instant in Mollweide view (a) near top of domain (O.95S0), (b) at mid 
depth (0.85Rq), and in (c) for an equatorial cut in longitude over full depth 
range. Strong plumes span the convection zone in the equatorial regions only 
within the nest of enhanced convection. The weaker cellular flows outside 
the nests are confined by shear to the upper reaches. 

motions is pronounced at all depths within the nest of active 
convection. 

The downfiowing plumes are influenced by the strong ra- 
dial shear and some break into multiple cells with radius even 
before the full blown nests of localized convection emerge, as 
is evident already in our simulation rotating at twice the so- 
lar rate (case G2). When the downflow networks only span a 
portion of the convection zone and experience a limited range 
of the full density stratification, the importance of compress- 
ible effects decreases. This has important consequences for 
the energetics of the simulations, particularly the radial ki- 
netic energy flux, as will be addressed in $6] In contrast, 
the downfiowing plumes in the polar regions experience much 
less shear from either the differential rotation or the relatively 
weak meridional circulations and continue to span the entire 
convection zone depth. 

3.4. Thermal Structuring 

In these rapidly rotating suns, the turbulent alignment of 
convection with the rotation axis leads to a net latitudinal 
transport of enthalpy, yielding a prominent latitudinal gradient 
of temperature. The resulting thermal structuring in case G5 is 
shown in Figure|6] presenting both the mean temperature pro- 
file and representative temperature fluctuations in a snapshot 
near the surface. In the latter, individual convection cells are 
associated with small fluctuations with amplitudes of a few 
K. Downflows are generally cool while upflows are relatively 
warmer. The enhanced enthalpy transport within the active 
nests of convection appears as positive temperature fluctua- 
tions in the equatorial region. 

Evident at high latitudes (Fig. |6j;) are broad spatial struc- 
tures (in addition to small-scale convection) which appear in 
the temperature fluctuations and are not readily visible in the 
maps of radial velocity (see Fig. at same instant). These 



structures are long lived and appear to be a separate phenom- 
ena from the nests of convection. The polar patterns propa- 
gate in a retrograde sense more rapidly than the differential 
rotation in which they are embedded, and though streaming 
wakes from the active nests print weakly into the polar re- 
gions, the polar patterns and nests appear to be distinct phe- 
nomena. The large-scale polar patterns are not evident in the 
slowly rotating cases (Gl and G2); in the most rapidly ro- 
tating cases this modulation attains a more complicated form 
than the two-lobed structure shown here. 

The zonally-averaged thermal structure (Fig. [6]2, Z?) is quite 
smooth and is characterized by warm poles and a cool equa- 
tor, with yet cooler mid-latitudes. In contrast, the mean en- 
tropy increases monotonically from equator to pole, due to 
effects of pressure. All of the more rapidly rotating cases 
have similar latitudinal thermal profiles, though the tempera- 
ture difference between equator and pole increases with more 
rapid rotation, as will be discussed further in §|4] In case G5, 
the latitudinal pole to equator temperature contrast is approx- 
imately 100 K throughout the convection zone. These latitu- 
dinal variations remain small at all rotation rates in compari- 
son to the spherically symmetric background T, which ranges 
from 2.7 x 10 5 K near the surface to 2.2 x 10 6 K near the bot- 
tom of the convection zone (as shown in Fig . [TJ . 

4. BUILDING A STRONG DIFFERENTIAL ROTATION 

4.1. Thermal Wind Balance 

Rapidly rotating systems are constrained by the Taylor- 
Proudman theorem to have minimal variations in flow dy- 
namics along the direction of the rotation axis. In strati- 
fied flows, gradients in density and pressu re contribute to 
baroclinic t erms in the vorticity equations (Ped loskylll982t 
IZahnlfT992|) which maintain flows that can break the Taylor- 
Proudman constraint. In our rapidly rotating suns, convec- 
tive plumes tilt towards the rotation axis as rotation effects in- 
crease. This results in latitudinal as well as radial transport of 
enthalpy and builds a latitudinal gradient of temperature and 
entropy. Such gradients arise naturally in a rotating convec- 
tive system even with uniform thermal boundary conditions. 
For a nearly adiabatic stratified, rotating, non-magnetized 
fluid it can be shown that in the limit of small Rossby num- 
ber and negligible viscous effects the zonal component of the 
vorticity equati ons reduces to the well known therma l wind 
balance (e.g.. lBrun & Toomrell2002tlMiesch et"aT1l2006l) : 



dz 2Q>rQ o <90' v 
where z is directed along the rotation axis and a hat denotes 
an average in longitude and time. We have further assumed 
that the turbulent pressure is negligible. 

From equation dD it is clear that departures from rotation 
constant on cylinders (as observed in the solar interior by he- 
lioseismology) can be maintained by a latitudinal gradient of 
entropy. The left and right hand sides of equation ([8]) are 
shown for case G5 in Figure [7] In the bulk of the convec- 
tion zone, the differential rotation profiles realized in these 
more rapidly rotating suns are substantially in thermal wind 
balance. Significant departures arise near the inner and outer 
boundaries (Fig. |7b) where Reynolds stresses and boundary 
conditions play a dominant role, as was found in earlier sim- 
ulations of solar convection ( Brun & Toomre 200l. 

Another striking property of the thermal wind balance is 
that increasing Qq leads to more cylindrical profiles of vj, un- 
less dS/d6 also adjusts with the rotation rate. In our more 
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FIG. 6. — Temperature structures within case G5. Mean latitudinal variations in temperature are shown relative to their spherical average T in (a) as contours 
with radius and latitude and (b) as cuts at fixed radii at the top (solid, Q.96Rq), middle (dashed, O.84K0) and bottom (dotted, 0.72Rq) of the domain, (c) 
Temperature fluctuations in a snapshot near top of domain (0.95Rq) relative to the mean structure in (a). 




FIG. 7. — Thermal wind balance achieved in case G5. (a) Gradients of vl 
along the rotation axis, dv^/dz, (b) the scaled latitudinal entropy gradient 
from the right-hand side of eq. {8), and (c) their difference, with contours in 
the latter rescaled to show the departures near the boundaries. The bulk of the 
convection zone is in thermal wind balance, but substantial departures arise 
near the top and bottom of the domain where Reynolds stresses dominate. 

rapidly rotating suns we find that the latitudinal gradients of 
temperature and entropy increase with more rapid rotation. 
The growth of AS (difference between the surface value of S 
at say 60° and the equator) with increasing rotation rate firj 
is shown in Figure [8] The latitudinal structure of entropy is 
always monotonic in these simulations, with lower entropy at 
the equator and higher entropy at the poles. Convection in 
these more rapidly rotating systems establishes stronger lati- 
tudinal gradients, but not enough in these simulations to main- 
tain the profiles unchanged. 

Accompanying the growth of AS is a growth in the latitu- 
dinal temperature contrast, as shown by the maximum tem- 
perature contrast in latitude near the stellar surface in Fig- 
ure [8] Typically, the maximal contrast occurs between the 
poles and latitudes of ±40°, as seen in Figure [6] for case G5 
with a contrast of about 100 K. In the rapidly rotating simula- 
tions, the primary flux balance in latitude is between thermal 
eddy diffusion KpT{dS/d9) and convective enthalpy trans- 
port C p p(v' g T'). Here convective transport moves warm ma- 
terial to the poles as the downflows align more strongly with 
the rotation axis while eddy diffusion works to erode the gra- 
dient. The meridional circulations appear to play a relatively 
minor role in maintaining the overall latitudinal entropy con- 
trast. 

4.2. Angular Momentum Redistribution 

In these simulations of stellar convection, complex cou- 
plings between rotation and convection build the profiles of 
differential rotation and meridional circulation. With stress- 




FIG. 8. — Increase with rotation rate in the latitudinal contrast of entropy, 
AS (plotted as diamonds) between equator and high latitudes at 0.96/f q . The 
more turbulent cases (G3a, G3b and G5b as labeled) have larger entropy con- 
trasts, in keeping with their generally stronger differential rotation. Blue tri- 
angles indicate the maximum temperature contrast in latitude at the upper 
boundary in each simulation. 

free boundary conditions at the top and bottom of the shell 
there are no net torques and thus the total angular momen- 
tum is conserved. Couplings between rotation and convection 
lead to a global-scale redistribution of angular momentum, 
resulting in the sustained flows of both differential rotation 
and meridional circulation. To assess the transport of angu- 
lar momentum in th ese systems we follow the approach of 
iMiesch et alj (|2008), examining the average radial and latitu- 
dinal angular momentum transport as detailed in their equa - 
tions (10)-(12) (see also lBrun & Toomrel2002tlMieschl 2005). 
The angular momentum fluxes from Reynolds stresses (RS), 
meridional circulations (MS) and viscous diffusion (VD) are 



F** = prsin( 



y r vj,)r+ i 



fMC . 



pc(ly r )?+(ve)t), 
pvr 2 sin 2 0Vfi, 



where 



C = rsin8 (il i) rsm9+ (v^)) 



(9) 
(10) 
(11) 

(12) 



is the specific angular momentum. 

The total radial and latitudinal fluxes of angular momentum 
are shown for case Gl and G5 in Figure [9] Here we have in- 
tegrated in co-latitude and radius respectively to deduce the 
net fluxes through she lls at various r adii and through cones at 
various latitudes (c .f . . lMieschll2005l) . The three major contri- 
butions arise from Reynolds stresses, meridional circulations 
and viscous terms. Velocity correlations lead to net angu- 
lar momentum transport by Reynolds stresses as convective 
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FIG. 9. — Time average of the integrated radial (F r ) and latitudinal (Fg) 
angular momentum flux for case Gl (a. b) and case G5 (c,d). Shown are the 
fluxes from Reynolds stresses (RS), meridional circulations (MC), viscous 
diffusion (VD) and their total. Transport by viscous diffusion remains com- 
parable in all cases, while the transport by Reynolds stresses and meridional 
circulations changes markedly with more rapid rotation. 

structures develop or ganized tilts and align partially with the 
axis of rotation (e.g.. iBrummell et al. 1998; Bru n & Toomrd 
l2002t iMieschet all 120081) . This alignment is particularly 
prominent in the fast downfiow lanes, and becomes stronger 
as rotation increases. 

Turning first to our solar case (Gl, Fig. [9jj, we see that 
in radius the meridional circulations and Reynolds stresses 
play similar and nearly equal roles in transporting angular 
momentum outward. The viscous flux meanwhile is negative 
and transports angular momentum inward, in keeping with 
the positive radial gradient of the differential rotation pro- 
file (eq. fTTT i. and the total flux in radius is nearly zero. The 
transport in latitude is somewhat different. Here meridional 
circulations combine with viscous fluxes to transport angu- 
lar momentum away from the equator and toward the poles 
(i.e., positive in the southern hemisphere and negative in the 
northern hemisphere). This tendency is opposed by Reynolds 
stresses, which continuously accelerate the equatorial regions 
and dynamically maintain the angular velocity contrast Af2 in 
latitude. 

The transport of angular momentum in our more rapidly ro- 
tating cases are all similar in form and are well represented by 
case G5 (Fig. |9j:, c/). In these more rapidly rotating stars, the 
radial balance is dominantly between the Reynolds stresses 
transporting angular momentum outward and viscous terms 
transporting it inward. The viscous transport is similar in 
magnitude to that of case Gl, though the radial boundary 
layers are now much narrower. The transport by Reynolds 
stresses is nearly twice as large, and this likely arises from the 
strong alignment of convective structures in both polar and 
equatorial regions. In these stars the weaker meridional circu- 
lations become relatively minor and disorganized players in 
the radial flux balance, moving angular momentum outward 
in some regions of the shell and inward in others. This oppos- 
ing behavior between the upper and lower convection zone 
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FIG. 10. — (a) Angular velocity contrast Af2 in latitude between equator 
and 60° (diamonds) and in radius across the shell at the equator (blue trian- 
gles). The more rapidly rotating cases appear to follow a power law, which 
for th e lat itudinal contrast is m = 0.3 and for the radial contrast is m = 0.4 (as 
in eg. I14t . (b) Relative latitudinal angular veloc ity c ontrast AC/f2eq, with 
the shown power law having n = —0.6 (as in eq. |13> . The scaling may vary 
with the path in parameter space, as suggested by cases G3a, G3b and G5b. 

arises from the meridional circulations breaking into multiple 
cells in radius. 

The balances achieved in the latitudinal transport in case G5 
(Fig. |9}i) are more complex. As the rotation rate has in- 
creased, the total viscous transport has remained nearly con- 
stant, with the significantly stronger gradients of angular ve- 
locity in the differential rotation profiles offset by the lower 
turbulent diffusivities dictated by our path through parameter 
space. That these two opposing actions should conspire to 
produce a nearly constant profile of viscous angular momen- 
tum transport is striking and not intuitive. This is particularly 
apparent when we examine the two other terms in the flux 
balance. The meridional circulations have reversed their role 
from our solar-like case Gl and now work with the Reynolds 
stresses to accelerate the equator and spin down the polar re- 
gions. The reduced contribution of the meridional circulations 
to the total balance arises as the flows become both weaker 
and multi-celled in radius and latitude. The smaller transport 
by Reynolds stresses appears to result from the destruction by 
radial shear of some of the downfiow plumes. 

4.3. Differential Rotation and Scaling with Rotation 

In analyzing our simulation results, it is the differential ro- 
tation established by the convection that may yield the most 
direct contact with observations. Stellar observations across 
the HR diagram indicate that differential rotation is a com- 
mon feature in many stars, particularly stars of spectral class 
F and later. In the sun, differential rotation has been mea- 
sured thro ughout the bulk of the convection zone (as re- 
viewed by iThompson et alJ 120031) . but at present for more 
distant stars only the surface differential rotation can be in- 
ferred. A variety of observational techniques have been em- 
ployed, ranging from photometric variability studies (e.g., 
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techniques (e. g.. iDonati et al.l 120031) a nd Fourier transform 
methods (e.g. jReiners & Schmittll2003l) . Typically, these ob- 
servations seek to measure the amount of angular velocity 
contrast at the stellar surface, denoted as Ail*, though what 
is being measured may be somewhat uncertain. Variations 
in AS1* have been found with both rotation rate and spectral 
type, however these quantities are correlated and in observa- 
tions to date it is diffic ult to disentangle their possible separate 
effects dReinersll2006l) . 

Different techniques measure fundamentally different trac- 
ers of surface differential rotation, either following variabil- 
ity of Ca emission (photometric), darkening from inferred 
starspot presence (photometric and Doppler imaging) or ro- 
tational broadening of absorption lines of unspotted stars 
(Fourier transform methods). Each technique also is most ap- 
plicable in only a limited region of stellar parameter space. 
As such, overlapping surveys are in short supply. Generally, 
most observations indicate that the relative shear, Af2»/f2*, 
depends on the stellar rotation rate SI* as a power law, though 
different surveys find different scalings for the differential ro- 
tation, expressed as 

AfL 

— iocO" (13) 



(e.g., n = -0.3 ± 0.1 inlDonahue et al.l ll996l and n = -0.34 ± 
0.26 in Reiners & Schmittl l2003T but n = -0.85 ± 0.10 in 
iBarnes et al.ll2.Q05b . Whereas some global models of convec- 
tion in more rapidly rotating stars have be e n conducted (e.g., 
RudieeretalJ [1998: iKuker&Stixl 12001 Ifuker & Riidiger 



2005MT these have been largely carried out in 2-D under 



the simplifying assumptions of mean-field theory. 

The amount of latitudinal shear observed at the surface Aft 
is an important quantity both for interpreting stellar observa- 
tions and for many dynamo theories. Here we define AJl 
more specifically as the difference in angular velocity be- 
tween the equator and say at 60° latitude, namely 



AS1 = Sl e q — ^60 OC SVq. 



(14) 



Going to higher latitudes yields comparable behavior. As 
shown in Figure [TUJ we find that Ail increases with rotation 
rate in our simulations, with m = 0.3 in the most rapidly ro- 
tating simulations. The radial shear also increases with more 
rapid rotation, and at the equator the difference between the 
mean angular velocity at top and bottom of the shell scales 
as m = 0.4 for the rapid rotators. Because m < 1, the relative 
shear AS7/il eq decreases with rotation rate f2o for the rapid 
rotators, in nearly a power law fashion for the path through 
parameter space explored here. The scaling exponent from 
Equation ( [T3l exhibited by these cases is n = -0.6, but this 
may be influenced by our choice in the scaling of diffusivi- 
ties with rotation. We are encouraged that our more turbu- 
lent cases G3b and G5b, with the same diffusivities at dif- 
ferent rotation rates, exhibit similar behavior. Our choice of 
low Prandtl num ber also has an effect on this scaling (see 
iBallot et"ai1l2007l) . Additionally, different treatments of the 
SGS unresolved flux, which has the most effect in the upper 
10% of the convection zone, can alter t he particular s caling 
law. The early simulations presented in iBrown et al.l (120041) 
and shown in Figure [2] which had a much thicker unresolved 
flux layer, had a scaling of n = -0.8 in the rapid rotation limit. 
We have found that normalizing Ail by f2rj rather than fi eq 
leads to a systematic offset for n of about -0.05 in the inferred 
scaling law. 
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FIG. 11. — Mean meridional circulations with latitude and radius for 
(a) case Gl, (b) case G5 and (c) case G10 with streamlines of mass flux 
overlaid. Colors indicate the sense (red counter-clockwise, blue clockwise) 

and magnitude of the meridional velocity (v„) = (v r )f + (vg)8. With more 
rapid rotation the meridional circulation cells align strongly with the rotation 
axis and weaken in amplitude, (d) Amplitude of the mean latitudinal com- 
ponent vg at the top of the simulation for case Gl (blue), G5 (black), and 
G10 (red), with regions of poleward and equatorward flow denoted. 

5. MERIDIONAL CIRCULATIONS AND SCALING WITH ROTATION 

The meridional circulations realized within our simulations 
are of significance since they can variously transport heat, an- 
gular momentum and even magnetic fields between the equa- 
tor and the poles, though the later are not included in the 
present simulations. Our time and longitude-averaged merid- 
ional circulation patterns are shown in FigureQT]for cases Gl, 
G5 and G10, depicted as streamlines of mass flux 

19* <9* 

rsm9(pv r )= — and rsm8(pv g ) = — (15) 

r oO or 

and averaged here over a period of at least 150 days. 

In our more rapidly rotating cases the meridional circu- 
lations have broken into several cells strongly aligned with 
the rotation axis (Fig. Q~TJ>,c), particularly in the equatorial 
regions. Weak connections between the equatorial and po- 
lar regions persist at the highest rotation rates studied, with 
organized flows along the tangent cylinder. These internal 
flows weaken with more rapid rotation. The meridional cir- 
culations are complex and time dependent, with large fluctua- 
tions around the statistically-steady states shown here, involv- 
ing variations comparable to or larger than the mean values 
themselves. The circulations are driven by small imbalances 
between relatively large forces and their nature is subtle. The 
variation of the meridional flows near the surface (0.96/?©) 
with rotation rate is shown in Figure QjJ/. The amplitude 
of the flows decreases substantially with more rapid rotation. 
Peak velocities drop from 22 ms _1 in case Gl to 14 ms"'in 
G5 and about 7 ms"'in G10. 

The total energy contained in these meridional circulations 
decreases quickly with more rapid rotation, as shown in Fig- 
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FIG. 12. — Volume-averaged kinetic energy density of the meridional cir- 
culations (MCKE) with rotation rate Qq. The MCKE is normalized by that 
energy in case Gl at the solar rate (2.5 X 10 4 ergs cm -3 ). The kinetic en- 
ergy of these circulations decreases with rotation rate; a power law scaling of 
Qq ' is shown for reference. 

ure Q~2] This drop in energy is independent of the detailed 
structure of the convection, showing no change in behavior at 
the transition to spatially modulated convection. In contrast 
to the energy contained in convection (CKE) and differential 
rotation (DRKE), the energy in the meridional circulations 
(MCKE) is much less sensitive to the level of turbulence in 
any particular simulation, as indicated by cases G3, G3a and 
G3b (detailed in Table |2). The meridional circulations remain 
important to the global-scale dynamics as their gradual redis- 
tribution of angular momentum contributes to the large angu- 
lar velocity gradients in latitude. Yet they are inefficient at 
transporting heat out of the star and at redistributing thermal 
material to maintain the latitudinal gradients of temperature 
and entropy (which correspond to the thermal-wind compo- 
nent of the achieved differential rotation). 

This finding is in striking contrast to the assumptions 
of many Babcock-Leighton dynamos, which often take the 
meridional velocity to scal e as v m oc £1 or v m oc logfi (e.g., 
ICharbonneau & Saarl 120011: iDikpati et al.1 1200 lb . In these 
models faster meridional circulations lead to shorter dy- 
namo cycles as surface flux is returned more rapidly to 
the tachocline. Currently, the observational data does not 
appear good enough to distinguish between the compet- 
ing models and we will have to await better measurements 
of the scaling between cycle period and rotation rate and 
possibl e observations of the meridional circulations them- 
selves (RempeJ l2008h . Existing 2-D mea n -field models 
of rapidly rotating suns ( Kuk er & Stixl l200lt iRiidiger et al.1 
1998; lKuker&Rudigerll2005ah also predict an increase of 
meridional circulation velocities with more rapid rotation. 
This is in contrast to their decrease in our simulations. 

6. ENERGY BALANCES AND FLUX TRANSPORT 

Convection is responsible for transporting the stellar flux 
emerging from the deep interior through the convection zone. 
In these simulations, the total luminosity L(r) and its compo- 
nents are 

F e +F k +F r +F u +F v =^=F t , (16) 

with 

F e = c p pv r T', F k =\pv r v 2 , (17) 

-rtf _-(iS =>| 

F r = -K r c p p^jr, F u = -K pT—, F v =-vV\ r , (18) 

where F e is the enthalpy transport by convective motions, F k 
is the kinetic energy flux, F r is the transport by radiation, F u 
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FIG. 13. — Variation of the energy fluxes with radius for cases Gl (dashed) 
and G5 (solid). Shown are the radiative luminosity (F r ), convective enthalpy 
transport (F e ), unresolved flux (F u ), kinetic energy transport (F k ), and the to- 
tal flux through the simulation (F,), all normalized by L,Q/(4irr 2 ). Case Gl 
has a large positive convective enthalpy flux throughout the domain, with the 
excess in luminosity largely balanced by an inward flux of kinetic energy. At 
higher rotation rates (as in G5), the kinetic energy flux is nearly zero through- 
out the domain. In both cases F v is negligible and not shown. 

is the unresolved SGS heat flux for parametrized transport by 
scales of motion below the resolution of our simulation and 
F v is the SGS viscous flux. Figure [13] shows the flux bal- 
ance with radius achieved in cases Gl and G5, averaged over 
horizontal surfaces and converted to relative luminosities. In 
the deepest layers the radiative flux becomes significant as the 
radiative conductivity steadily increases with depth. By con- 
struction it suffices to carry all of the imposed flux through 
the lower boundary, where the radial velocities and convec- 
tive flux vanish. A similar role is played near the top of the 
convection zone by the sub-grid scale transport which yields 
F u . The main role of F u is to transport energy outward through 
the impenetrable upper boundary where the convective fluxes 
vanish and the remaining fluxes are small, thereby avoiding 
the building of strong superadiabatic radial gradients there. 

The functional form of Ko(r) is chosen so that the entire 
stellar luminosity will be transported at the surface of the 
simulation by_,F„. A subtlety of this treatment of the SGS 
flux lies in dS/dr. In the more rapidly rotating simulations, 
we find that convection is less able to establish an adiabatic 
profile throughout the convection zone. Instead, much of the 
convection zone remains slightly superadiabatic (dS/dr < 0). 
This property is also realized in local-domain simulations of 
rapidly rotating convection, where more rapid rotation leads 
to enhanced horizontal mixing through vortex interactions 
and a resulting decrease in enthalpy transport as vertical ve- 
locities and therma l perturbations become more decorrelated 
(lJulien et alJ Il996t iBrummell etail 119961) . The change of 
dS/dr with rotation rate is shown in Figure [TJ; for four of 
our simulations. Fixing the amplitude and structure of Ko(r) 
across simulations leads F u to influence a greater portion of 
the convection zone at more rapid rotation, as indicated by 
the slight growth of F u in case G5. Though this effect be- 
comes stronger as our simulations rotate more rapidly, at no 
point in these simulations does F u transport more than 10% of 
the total luminosity at mid convection zone. 

In all these simulations, the strong correlations between ra- 
dial velocities and temperature fluctuations yield the enthalpy 
flux F e , which dominates the energy transport at mid convec- 
tion zone. Both warm upflows and cool downflows serve to 
transport flux out of the star, and the two carry comparable 
amounts of flux to one another in the rapidly rotating sim- 
ulations. In going to the more rapidly rotating simulations, 
we find that the average convective enthalpy flux through the 



Rapid Rotation and Modulated Convection 



13 



TABLE 2 
Flux Balances and Energies 



Case 


J? a 


F a 






CKE b 


DRKE b 


MCKE b 


AT c 
J max 


Gl 


1.014 


1.140 


-0.118 


-0.190 


3.28 


2.26 


0.025 


5.50 


G2 


1.300 


0.684 


-0.088 


0.046 


2.64 


13.2 


0.015 


28.0 


G3 


1.349 


0.628 


-0.077 


0.071 


2.40 


20.5 


0.011 


53.5 


G4 


1.327 


0.631 


-0.065 


0.073 


2.21 


25.5 


0.009 


78.5 


G5 


1.329 


0.625 


-0.079 


0.078 


2.11 


30.1 


0.007 


107 


G7 


1.298 


0.581 


-0.097 


0.080 


1.69 


38.7 


0.005 


171 


G10 


1.236 


0.623 


-0.093 


0.090 


1.51 


47.1 


0.003 


271 


G3a 


1.268 


0.655 


-0.071 


0.068 


2.73 


27.7 


0.012 


58.7 


G3b 


1.101 


0.780 


-0.065 


0.043 


3.34 


37.9 


0.013 


62.4 


G5b 


1.172 


0.668 


-0.047 


0.081 


2.44 


57.4 


0.008 


134 



a Average convective enthalpy (F e ) and kinetic energy (Ft) fluxes at mid-layer scaled 
by the solar flux, shown for polar (latitudes above ±60°) and equatorial (from ±30°) 
regions. b Kinetic energy density relative to the rotating coordinate system, for convection 
(CKE), differential rotation (DRKE) and meridional circulations (MCKE), averaged over 
the full shell and over ~ 300 days; units are 10 6 erg cirT 3 . c Maximum temperature contrast 
at 0.96Rq in K, typically occurring between pole and ±40° 



polar regions is greater than that through the equator. This 
is shown in Table [2] In case Gl (at the solar rate) slightly 
more flux is transported through the equator than the poles, 
but as the rotation rate increases significantly more flux is 
transported through the poles than the equator. This latitudi- 
nal variation becomes somewhat weaker in our more turbulent 
cases (G3a, G3b and G5b), but remains present. Interestingly, 
cases G7 and G10 follow similar trends in their equatorial en- 
thalpy transport, despite the emergence of strong nests of con- 
vection at the equator and the suppression of convection in the 
rest of that region. 

In case Gl, convection transports slightly more than the so- 
lar luminosity. This over luminosity is balanced by an inward 
transport of kinetic energy, which is primarily due to com- 
pressible effects and the transport of v 2 r within strong down- 
flows that sp an the convection zo ne and feel the full density 
stratification dHurlburt et alll9 86). In the more rapidly rotat- 
ing cases, the prominent differential rotation shears apart the 
convection cells. The downflows have lost much of their co- 
herence and only the strongest downflows within the nests of 
localized convection survive to span the full convection zone. 
Individual downflows thus feel less density stratification and 
compressible effects become less important, leading to a bal- 
ance in the transport of v 2 . between the upflows and down- 
flows. Instead, as shown in Table[2]the sense of Fj reverses in 
the equatorial regions, becoming dominated by the outward 
transport of v\. The polar regions remain largely unchanged 
and the total across spherical surfaces is nearly zero. 

Volume-averaged energy densities for a selection of our 
simulations are shown in Table [2] At the solar rotation 
rate, convective kinetic energy (with kernel jpv 1 and labeled 
CKE) and the kinetic energy in the average differential ro- 
tation (ip(v0) 2 , DRKE) are comparable. As the rotation 
rate is increased, DRKE grows strongly and convective en- 
ergy decreases slightly, leading DRKE to dominate the total 
energy budget. This is true even in our significantly more 
turbulent solutions. The energy in meridional circulations 
(|p((v r } 2 + (vg) 2 ), MCKE) is always small, and decreases in 
both magnitude and percentage of the total energy with more 
rapid rotation. 

7. ACTIVE NESTS OF CONVECTION 

The emergence of localized nests of convection at higher 
rotation rates is a striking feature that calls out for an expla- 



nation. In many ways, it is quite surprising that convection 
chooses to be confined to narrow intervals in longitude, but 
such states have also been realized in a number of other dy- 
namical systems. Generally the appearance of nests is a chal- 
lenge to explain in detail, yet the onset of spatially modulated 
states which may be their precursor is better understood. 

7.1. Spatially Localized Convection in Other Settings 

The phenomena of spatially localized convection has a rich 
history, variously appearing in laboratory experiments and 
numerical simulations. Much interest in confined states of 
convection was motivated by t he discovery of such states in 
binary fluid convection (e.g.. lAnderson & Behringer 119901: 
Kolod ner & Glazier! 119901: iNiemela et alJ 119901: ISurko et all 
1199 lb . where traveling waves of convection appear via 
subcritical Hopf bifurcations and near onset are seen to 
evolve into traveling patches of convection separated by re- 
gions of nearly quiescent fluid. From a theoretical per- 
spective, these confined states near onset are accessible 
to weakly nonlinear theory and considerable progress has 
been made in unde rstanding their nature (e.g ., Rieck el 11992b 
i Barten et all 119951: IBatiste & Knobiochl 120051; IBatiste et alJ 
l2006tlBurke & Knob loch 2007). The confined states found in 
binary convection differ from our active nests of convection in 
several respects. The most important is that in binary convec- 
tion there is no net vertical transport of solute. The confined 
states instead pump solute horizontally and create regions of 
stable vertical stratification in the quiescent regions. A possi- 
bly related phenomena is that of localized sta tes in magneto- 
convection, as studied by Blanchfl owerl d 1999b and in 3-D by 
iBlanchflower & Weiss! (120021) . Here single convective cells 
(called "convectons") formed in a region of initially uniform 
strong vertical magnetic field by a process of flux expulsion. 
Convection within the localized states was strong and was en- 
tirely suppressed in surrounding medium. These convectons 
could contain several convective cells and were generally sta- 
tionary, though some solutions exhibited time dependent be- 
havior. Recently, progress has bee n made addre ssing these 
systems in approximate 2-D models (lDawesll2007l) . 

Confined states are also realized in other doubly diffusive 
systems, as in theoretical studies of thermosolutal convection 
( Knob iocrTeTaT] [1985 iDeane et alJ [19871 fl988l ISpiiia et alJ 
119981) . In the latter studies a variety of traveling convec- 
tive patches were found, and in these the convective trans- 
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port of heat and solute was enhanced compared to that in 
the interpatch regions. In all cases the patches propagated in 
the same direction as the individual convective cells, though 
more slowly. Such behavior persisted for long periods of time. 
These localized states occurred well above the onset of con- 
vection, and convection continued in the interpatch regions. 
There may also be analogues in convection within the Earth's 
atmosphere, where deep convection in the tropics tends to be 
organized on global scales into regions of locally enhanced 
convection which propagate in a prograde sense. These orga- 
nized convective structures are called the Madden-Julien Os- 
cillation and appear to have their origin in the coupling of 
convective motions with equatorially tr apped wave s (perhaps 
Rossby or Kelvin waves; see review by Zhang 2005). 

The spatially modulated states in our simulations of stel- 
lar convection exist at Rayleigh numbers far above the onset 
of convection. Spatially modulated states in this parameter 
regime have also been observed in geophysically motivated 
3-D Boussinesq simulati ons of convection within a thick , 
rotating spherical shel l (iGrote & Bussd 120001: iBussd 120021: 
iBusse & Simitevl[2005l) . In these studies, spatially localized 
states emerged at moderate Rayleigh numbers, involving an 
equilibrium between the shearing flow of differential rotation 
destroying convective eddies and the Reynolds stresses in the 
convection driving the differential rotation. These effects led 
to localized states where convection occupied a limited por- 
tion of the domain and the region outside the convective patch 
was filled with quiescent streaming flow and almost no radial 
motions. In the quiescent regions the thermal gradients be- 
come increasingly unstable until they are advected back into 
the patch where they help sustain the convective eddies. The 
patches in these geophysical simulations moved slowly ret- 
rograde and persisted up to Rayleigh numbers of about 10 6 . 
Beyond this point the differential rotation became so strong 
that no sustained convection was possible. Instead the system 
began to behave as a relaxation oscillator, with short bursts 
of convection temporarily building a strong differential rota- 
tion which then sheared out the convective eddies. Convec- 
tive transport remained suppressed until the shear of differ- 
ential rotation decayed viscously, whereupon a new burst of 
convection would begin the process anew. In all cases with 
localized convection, significant oscillations were seen in the 
kinetic energies of bo th differential rotation and convection 
(IGrote & Bussdl2000l) . 

Similar states have also been realized in anelastic sim- 
ulations of stellar convection in spherical shells for a ro- 
tating younger sun with a much thicker convection zone 
(iBallot et al.l 120061 120071) . The spatially modulated states 
found there appear in the equatorial convection for models 
with low Prandtl number (Pr = 0.25, as in the models of this 
paper). Localized states turn into bursty, vacillating convec- 
tion at large Taylor n umbers (Ta > 10 9 ), much like those in 
IGrote & Bussd d2000l) . Localized states observed in thick 
convective shells (with typical aspect ratios x = '"bot/'top = 
0.58) differ in many important respects from the states we find 
in our relatively thin shells of convection (x = 0.76), most no- 
tably in their temporal stability, which we will next discuss. 

7.2. Properties of the Active Nests 

Our active nests of convection appear first as regions of 
mildly enhanced convective amplitude. As the rotation rate 
increases, convection in the equatorial regions gradually be- 
comes more localized (as in Fig. [3]) and eventually is present 
only within the active nest. In some of our systems we ob- 



serve two nests or patches in longitude (such as case G5) and 
in some only a single nest (as in case G10). Generally, con- 
vection at the highest rotation rates possesses only a single 
nest, and for moderate rotation rates the system can alternate 
between two-nest states and single-nest states, suggesting that 
several attractors exist within the phase space. 

To study the temporal evolution of our convective patterns 
in more detail, we employ time-longitude maps as shown in 
Figure [14] Here radial velocities are sampled at the equator 
for all longitudes, considering case G5 at three depths (near 
the top, middle and bottom of the convection zone) and doing 
so over an interval of 260 days (or about 45 rotation periods). 
In these maps, structures propagating in a prograde fashion 
relative to the frame of reference are tilted to the upper-right 
and patterns propagating in a retrograde sense tilt to the upper- 
left. To construct these maps we have chosen a frame of ref- 
erence propagating in a prograde sense relative to the bulk ro- 
tation frame (Qq = 1 -3 x 10~ 5 rad s" 1 or 2070 nHz), with con- 
stant relative angular velocity 6.75 x 10~ 7 rad s -1 (107 nHz, 
for a total rotation rate of 1.052 ilo). This corresponds to 
the propagation rate of the modulated convection pattern, and 
thus the nests appear stationary in this frame. In the f^o refer- 
ence frame the nest takes about 108 days to complete one lap 
around the equator, thus for the time interval shown in Fig- 
ure[l4]the nests travel about 870° in longitude and completes 
about 2.4 circuits around the equator. The active nests of con- 
vection can persist for extremely long periods of time. The 
two nests visible in Figure [14] remain coherent for over 5000 
days of simulated time. 

Individual convective upflows and downflows appear as 
streaked red and blue regions respectively. In the upper con- 
vection zone (0.95Rq), the convective cells propagate more 
rapidly than the nests of convection. Here they overtake the 
nests from behind (from smaller longitudes) and then slowly 
swim through at a reduced speed. Within the nests convec- 
tive cells collide and interact, and the strongest survive to 
emerge from the front of the nest, where they speed up as 
they then propagate through the more quiescent region. When 
this occurs, typically a small wave train comprised of 2 to 3 
upflow/downflow pairs escapes, and as they enter the quies- 
cent regions these convective cells speed up to once again out- 
pace the zonal flow. In the lower convection zone {Q.73Rq), 
the nests propagate more rapidly than the convective cells, 
and the individual upflows and downflows appear as strong 
retrograde-directed streaks. Also visible in the lower convec- 
tion zone are low-amplitude velocity structures of rapid retro- 
grade propagation. They appear to be the weak equatorward 
extension of the large-scale (retrograde rotating) polar pat- 
terns evident in Figure[6] At mid-convection zone (0.85/?©), 
nearly all vertical flow is occurring within the nests of active 
convection, though the strongest cells outside the nest in the 
upper convection zone manage to weakly print down to this 
intermediate depth. 

The typical angular velocities of the differential rotation, 
the individual convective cells and the active nests of convec- 
tion are shown in Figure [l5]for case G5 and detailed for sev- 
eral cases in Table [3] The nests of convection propagate at a 
constant angular velocity at all depths in the convection zone 
and over the entire range of latitudes (±30°) where they are 
present. In contrast, the angular velocity fl of the differential 
rotation varies substantially with radius and latitude. At all 
depths, the individual convective cells propagate more rapidly 
than the zonal flow of differential rotation which they drive. 
Because the nests of convection propagate at an intermediate 
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FIG. 14. — Time-longitude map for case G5 of radial velocity v, sampled at three depths (a) near top, (b) middle and (c) bottom of layer. These samples 
are extracted at the equator, using a reference frame tracking the nests and starting from a mature state in the simulation. Two nests are visible, and individual 
convective cells appear as paired upfiows and downflows, which propagate slightly faster than the mean zonal flow that they establish and thus pass through the 
nests, (d, e, f) Time averages with longitude of entropy fluctuations S and convective kinetic energy density K in these samples. 
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TABLE 3 

Angular Velocities of Various Structures 



Case 



n 



^eq.bot 



Oct 



Gl 


2.60 




0.341 


0.060 


0.535 


0.180 


G3 


7.80 


0.511 


1.094 


0.220 


1.176 


0.330 


G5 


13.00 


0.675 


1.319 


0.197 


1.381 


0.288 


G10 


26.00 


0.830 


1.497 


0.140 


1.623 


0.280 


G3a 


7.80 


0.690 


1.295 


0.265 


1.421 


0.373 


G3b 


13.00 


0.670 


1.555 


0.250 


1.590 


0.413 


G5b 


13.00 


0.750 


1.778 


0.246 


1.800 


0.430 
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FIG. 15. — Various angular velocities with radius for features in case G5. 
The mean shearing zonal flow of differential rotation Q is denoted for all 
depths at the equator (solid curve) and at ±15° latitude (dotted). The propa- 
gation rate of the active nests (diamonds) is constant across the entire latitude 
and radial range where they appear (±30°). The average velocity of indi- 
vidual convective cells at the equator (asterisks) is more rapid than the mean 
zonal flow they establish. Their propagation is faster than the nests near the 
top and slower near the bottom. The global rotation rate f2o is 2070 nHz 
(marked). 

prograde rate, they experience a head-wind from the differ- 
ential rotation in the deep convection zone and a tail-wind 
near the surface. Despite this strong radial shear, the nests re- 
main coherent across the entire convection zone for long peri- 



NOTE. — All angular velocities in /^rad s and, except for 
the frame rate f!o, are given relative to Slo. The differential rota- 
tion at the equator f2 K ] is measured at 0.95Kq (top) and 0.73Rq 
(hot). The mean propagation rate of individual convective struc- 
tures Q c is measured at the same depths in time-longitude maps 
of v r and has a typical variance of ±0.04 ^trad s -1 . 

ods compared with either the lifetime of individual convective 
cells (10-30 days), the rotation period of the star (5 days) or 
typical thermal and viscous diffusion times (r K = 910 days, 
r„ = 3640 days, both at mid-depth). The time for the differ- 
ential rotation to lap the nests at the equator is about 112 days 
near the surface and 143 days at the bottom of the shell. 

At the higher rotation rates, single stable nests of convec- 
tion dominate the equatorial region. Time-longitude maps 
are shown in Figure [16] for the equatorial radial velocity at 
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FIG. 16. — As in Fig. [14] time-longitude map for case G10 of radial velocity 
(a) near top and (b) near bottom of layer. Here a single nest is realized, with 
individual convective cells almost entirely confined to the nest, though they 
continue to propagate within it. (c, d) Time averages 5 and K constructed in 
the same reference frame. 

two depths in case G10, our most rapidly rotating simula- 
tion. Within the nest, convection remains vigorous, with the 
strongest downflow networks still spanning the entire depth 
of the convection zone. This nest again propagates in a pro- 
grade sense relative to fio, with constant relative angular ve- 
locity 8.3 x 10~ 7 rad s" 1 (132 nHz, for a total rotation rate of 
1 .032 flo). Over the interval shown the nest completes almost 
3 circuits of the equator relative to the Oq reference frame. 
Individual convective cells continue to swim through the nest, 
moving more rapidly near the surface and more slowly in the 
deep convection zone. The very strong radial shear prevents 
all but the strongest of downflows from spanning the full con- 
vection zone. In the region outside the nest, convection is 
strongly suppressed and the main features are the weak flows 
associated with the retrograde propagating polar pattern. In 
the upper convection zone, occasional weak convective dis- 
turbances appear upstream of the nest. As these fluctuations 
enter the nest they grow in amplitude. 

Our nests of active convection may owe their existence to a 
competition between the shearing action of differential rota- 
tion acting on the individual convective eddies and Reynolds 
stresses within the convection helping to mai ntain the strong 
zonal flow s. Unlike the syste ms studied by iGrote & Bussd 
(2000) and Ball ot et ail d2007t) . our patchy convection is not 
accompanied by relaxation oscillation behavior or large ex- 
changes between the kinetic energy in convection and in the 
differential rotation. Rather, our nests are not bursty in time 
and instead persist for long periods in quasi-steady states. 
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FI G. 1 7. — Snapshot of nest structure in case G10, at relative day 221 of 
Fig. 1161 (a) Radial velocities in an equatorial cut shown for all longitudes 
and radii, (b) Companion entropy fluctuations S about their spherical means 
S (in cgs units), with high entropy fluid in red tones, (c) Radial velocity in 
the equatorial region shown near the top of layer. 

This is true even for our most rapidly rotating case consid- 
ered here (G10) and at much higher turbulent driving (G5b), 
though these simulatio ns exist at Taylor n umbers below the 
threshold suggested bv lBallot et ail d2007l) (Ta > 10 9 ). Cou- 
pling between equatorially trapped waves and convection may 
also have a role, but the contribution of this coupling to spatial 
localization has been difficult to elucidate. 

7.3. Detailed Structure of an Active Nest of Convection 

We focus here on the structure of the nests so evident in 
Figures FBI and [TSI for cases G5 and G10, but devote par- 
ticular attention to the single nest realized in the latter as a 
representative case. The active nests extend throughout the 
depth of the convection zone. This is illustrated most clearly 
in Figure fTTl;. showing the vertical profile of radial velocities 
with longitude in a cut at the equator. Convection is broken 
into multiple cells, one set above 0.9R Q and another below 
O.87?0. Only within the nest do strong downflows span the 
convection zone. The nest is embedded in a region of strong 
latitudinal and radial zonal shear (as shown by the tilted na- 
ture of its structure), yet the pattern propagates at a constant 
angular velocity at all depths and latitudes. At the surface 
and near the equator the nest propagates more slowly than 
the zonal wind, whereas at the base of the convection zone 
it propagates more rapidly. Individual convective structures 
swim still more rapidly and both enter and exit the region of 
enhanced convection. As such, the nest experiences a strong 
retrograde flow at the base and a strong prograde flow near the 
surface. 

The thermal structure of the nests is revealed by their en- 
tropy fluctuations, as shown in Figure [l7t>. In the upper con- 
vection zone, the mean zonal flow advects low entropy fluid 
into the nests from the left side. This fluid is then swept 
away by intermittent downflows and replaced by higher en- 
tropy fluid from below. In the lower convection zone, higher 
entropy fluid is swept into the nest from the right and lower 
entropy exits to the left. At mid convection zone, regions out- 
side the nest remain convectively unstable, but only weak ra- 
dial motions are driven here. 

The mean longitudinal structure of the nests can be assessed 
by forming temporal averages of various properties in a frame 
co-rotating with the nests of convection (these averages de- 
noted by a tilde). This is done at the equator for entropy 
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FIG. 18. — Mean circulations associated with the nest in case G10. (a) 
Time averaged radial (colors) and horizontal flows (arrows) near the top and 
(b) near the bottom of the shell, in a reference frame tracking the nest. The 
strong zonal flow of differential rotation has been subtracted at each latitude 
and is shown by the solid curves (both scaled by 125 m s~'arrow length, with 
zero velocity relative to the nest indicated by the vertical line). Within the nest 
strong eddy currents partially decelerate the flow, while outside the streaming 
mean zonal flow dominates. 

S and convective kinetic energy density K (with same form 
as CKE) for case G5 in Figure [PR/- f and for case G10 in 
Figure [16b , d. There are distinctive differences between the 
leading (to the right) and trailing (to the left) portions of the 
profiles, with S showing a gradual rise and steeper drop in 
going to decreasing longitudes (from right to left). The enve- 
lope of S is largely similar in form at the top and bottom of 
the convection zone. In contrast K, which traces the fluctuat- 
ing velocities of convection, changes form with depth, being 
skewed in the direction of the mean zonal flow. At the top of 
the convection zone this sense of skew is toward the right, and 
at the bottom it is toward the left. 

These nests must have gradual circulations associated with 
them, though these are a challenge to discern as they are much 
weaker than the strong mean zonal flows. We can get some 
sense of them by averaging the vertical and horizontal flows 
in the surroundings of a nest while tracking it over long time 
intervals. Shown in Figure [18] are the slow circulations asso- 
ciated with the nest in case G10 near the top and bottom of 
the convection zone. These flows have been averaged over 
a period of 615 days starting at the beginning of the interval 
shown in Figure [16] The flows suggest a systematic zonal 
circulation from fore to aft in the upper convection zone and 
toward higher latitudes and diverging from the nest at depth, 
though to achieve this we have subtracted the fast mean zonal 
flow which varies with latitude as shown. The mean upflows 
and downflows in the nest are quite weak, with amplitudes 
of a few m s" 1 , as compared with the convection which can 
have peak velocities of a few hundred ms" 1 . Test particles 
released in the flow would not simply circulate according to 
these mean circulations and would be instead swept along by 
the strong mean zonal flow and of course by the vigorous con- 



vection cells that propagate through the nest. Some test par- 
ticles would encounter strong downflows and would be swept 
down through the nest to the bottom of the convection zone 
while others would likely be carried horizontally out the nest 
and remain at a similar depth. Yet Figure [18] indicates that a 
weak large-scale circulation is realized, and this may serve to 
slowly pump fluid through the system. These mean circula- 
tory flows may also serve to inform analytic models of such 
nests of convection. 

8. CONCLUSION 

When stars like our sun are young they rotate much more 
rapidly than the present sun. In these stars rotation must 
strongly influence the convective motions and may lead to 
stronger global-scale dynamo action. We have explored here 
the effects of more rapid rotation on global-scale convection 
in simulations of stars like our sun. The mean zonal flows of 
differential rotation become much stronger with more rapid 
rotation, scaling as AS1 oc ilj} 3 or as AJ7/S! eq oc ^q 0,6 - In 
striking contrast, the meridional circulations become much 
weaker with more rapid rotation, and the energy contained in 
them drops approximately as ^q ' 9 . Accompanying the grow- 
ing differential rotation is a significant latitudinal tempera- 
ture contrast, with amplitudes of 100 K or higher in the most 
rapidly rotating cases. The maximum temperature contrast 
near the surface occurs between the hot poles and the cool 
mid latitudes at about ±40°. If this latitudinal temperature 
contrast prints through the vigorous convection at the stellar 
surface, it may appear as an observable latitudinal variation 
in intensity. The thermal contrast would presumably persist 
for long periods compared to stellar activity, offering a way to 
disentangle this intensity signature from that caused by spots 
of magnetism at the stellar poles. 

These simulations are entirely hydrodynamic and this pro- 
vides a major caveat to our findings on the scaling of differ- 
ential rotation and latitudinal temperature contrast with rota- 
tion rate Oo- Prior MHD simulations of stellar convection 
have demonstrated that in some parameter regimes strong 
dynamo-generated magnetic fields can react back strongly 
on the differential rotation, acting to lessen a ngular ve- 
locity contrasts or largely ehm inate them (e.g. , iBrun et all 
2005: Feather stone et alj 120071: iBrowningl 120081) . It is un- 
clear whether the scaling trends identified here for differ- 
ential rotation as a function of f2o will survive in the pres- 
ence of dynamo action and magnetic fields. Likewise, mag- 
netic fields may lessen the strong temperature contrasts real- 
ized here, doing so through their feedback on the convective 
flows and energy transport. We expect that the weaker merid- 
ional circulations may be less affected by magnetic feedbacks, 
and thus the prediction that meridional circulations lessen in 
energy and amplitude with more rapid rotation may be of 
greater significance though harder to confirm observationally. 
Weaker meridional circulations in more rapidly rotating stars 
will have a strong impact on many theories of stellar dynamo 
action, including the Babcock-Leighton flux-transport model 
favored for solar- type stars. We have begun simulations to ex- 
plore the dynamo action realized at various rotation rates, and 
its impact upon the flows des cribed here. Preliminary results 
appear in iBrown et alJ (120071) and more detailed results will 
be forthcoming shortly. 

A striking feature of these simulations is the emergence 
of a pattern of strongly modulated convection in the equato- 
rial regions. These nests of active convection are regions of 
enhanced convective vigor and transport which propagate at 
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rates distinct from either the mean zonal flows of differential 
rotation or the individual convective cells. In the most rapidly 
rotating systems, such as case G10, convection at the equa- 
tor is entirely dominated by motions inside the nest with only 
very weak radial motions present in the regions outside the 
nest. Though their impact on the convection is most obvious 
in the rapidly rotating limit, we find some evidence for weak 
modulation even in our more slowly rotating cases. 

All of our simulations stop short of the turbulent stellar sur- 
face, and it is thus difficult to estimate how these nests of ac- 
tive convection may affect stellar observations in detail. The 
convective velocities associated with the nests are small com- 
pared to the nearly supersonic flows in stellar granulation, 
and in the sun such global-scale convective structures have 
evaded direct detection despite intensive searches throughout 
the near-surface layers by helioseismology. The extremely 
localized nests found in our most rapidly rotating cases may 
however influence the thermal stratification and thus convec- 
tive vigor in the near surface regions, as most of the flux at the 
equator is transported through a narrow range of longitudes. 
These nests may act as traveling hot spots with enhanced con- 
vection even in the surface layers where the higher emerging 
flux escapes the system. 

These spatially localized states of convection may also have 
some bearing on the active longitudes of magnetic activity ob- 
served in the sun, if the enhanced pummeling of the tachocline 
by the convection within the nest preferentially destabilizes 
magnetic structures within the tachocline that then rise to the 
surface. Initial dynamo simulations indicate that nests of con- 
vection can coexist with magnetism in portions of parame- 



ter space. Thus their strongest signature is likely to emerge 
in magnetic stars, where magnetic fields threading the bulk 
of the convection zone may be concentrated in the nests and 
mimic giant, propagating starspots which survive for very 
long epochs in time. If the nests lead to active longitudes of 
enhanced magnetic activity in rapidly rotating stars, we might 
expect these long-lived magnetic structures to propagate at a 
rate different from the stellar rotation rate as measured at the 
surface or from the stellar differential rotation. 

We recognize that our simulations remain separated by 
many orders of magnitude from the parameter space of real 
stellar convection. As such, we must be cautious with our 
interpretations of the overall dynamics. However, we have 
found that these nests of convection are a robust feature over 
a range of parameters and that they are able to persist as enti- 
ties for as long as we could pursue their modelling. Thus one 
should be prepared to consider the possibility of their pres- 
ence also in real stellar convection zones, where they may 
appear as long-lived propagating features. 
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